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The present study is to develop a new user-defined function using artificial neural networks intent Com- 
putational Fluid Dynamics (CFD) simulation for the prediction of water-vapor multiphase flows through fuel 
assemblies of nuclear reactor. Indeed, the provision of accurate material data especially for water and steam 
over a wider range of temperatures and pressures is an essential requirement for conducting CFD simulations 
in nuclear engineering thermal hydraulics. Contrary to the commercial CFD solver ANSYS-CFX, where the 
industrial standard IAPWS-IF97 (International Association for the Properties of Water and Steam-Industrial 
Formulation 1997) is implemented in the ANS YS-CFX internal material database, the solver ANS YS-FLUENT 
provides only the possibility to use equation of state (EOS), like ideal gas law, Redlich-Kwong EOS and piece- 
wise polynomial interpolations. For that purpose, new approach is used to implement the thermophysical prop- 
erties of water and steam for subcooled water in CFD solver ANSYS-FLUENT. The technique is based on 
artificial neural networks of multi-layer type to accurately predict 10 thermodynamic and transport properties 
of the density, specific heat, dynamic viscosity, thermal conductivity and speed of sound on saturated liquid and 
saturated vapor. Temperature is used as single input parameter, the maximum absolute error predicted by the 
artificial neural networks ANNs, was around 3%. Thus, the numerical investigation under CFD solver ANSYS- 
FLUENT becomes competitive with other CFD codes of which ANSYS-CFX in this area. In fact, the coupling 
of the Rensselaer Polytechnical Institute (RPT) wall boiling model and the developed Neural-UDF (User Defined 
Function) was found to be useful in predicting the vapor volume fraction in subcooled boiling flow. 
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I. INTRODUCTION 


The present study is concerned with artificial neural net- 
works support to the ongoing study on nucleate boiling. 
The focus is placed on numerical simulations of the pre- 
conditioning of two phase flows in a circular tube with heated 
wall under pressurized conditions. Computational fluid dy- 
namic techniques show increasing promise for the simulation 
of subcooled nucleate boiling [1]. The commercial Com- 
putational Fluid Dynamics (CFD) solver ANS YS-FLUENT 
is employed as the computational platform and the Eulerian 
method in combination with RPI boiling model is used, where 
the phases are assumed to be interpenetrating continua. This 
model is the most general and complex of all multiphase flow 
models. The advantage of ANSYS-FLUENT is the number 
of turbulent model adapted to special case. On the other hand, 
CFX was the go-to code for turbomachinery. 

Contributions are made towards an improved flow boil- 
ing model under Pressurized Water Reactor (PWR) condi- 
tions by carrying out detailed CFD analyses based on ex- 
perimental data of void distributions in heated tubes. When 
a heated surface exceeds the saturation temperature of sur- 
rounding coolant, boiling on the surface becomes possible. 
Bubbles formed on the heated surface depart the surface and 
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are transported by the bulk fluid, such that a condition of two- 
phase flow is said to exist. Depending on degree of subcool- 
ing and length of the heated tube, the bubbles may or may 
not condense and collapse prior to exiting the tube. In sub- 
cooled boiling this process results in further heating of the 
fluid toward the saturation temperature. In saturated or bulk 
boiling, bubbles can be transported along the entire length of 
the heated tube without collapsing [2]. 

In Rensselaer Polytechnical Institute (RPI) wall boiling 
model or heat partitioning model, implemented in the solver 
ANSYS-FLUENT version 14.5.0, the overall heat flux from 
the heated wall to the two-phase flow includes convective, 
quenching and evaporation heat flux. Furthermore, the heat 
flux partitioning model associates each of the heat flux con- 
tributions with a dimensionless wall area ratio, so as to define 
the ratio between heat flux contributions [1]. An important 
aspect is the knowledge and understanding of the fluid behav- 
ior in process of heat transfer. The thermophysical properties 
of water and steam are indicated as the fundamental proper- 
ties that govern the process [3]. 

There are several ways to focus on the thermodynamic and 
transport properties: through different state equations and by 
a number of approximation functions. In the early 1990s, 
the appearance of safety analysis codes of nuclear facilities 
boosted the search on accuracy and expanded the range of 
water and steam properties. Then, the IAPWS new standards 
leveraged data and technological advances towards new for- 
mulations of thermodynamic properties [4]. 
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Artificial neural networks are being extensively used for 
predicting the 10 thermodynamic and transport properties 
mentioned above. Artificial neural network (ANN) is an ad- 
vanced mathematical tool which determines the network out- 
put based on available experimental information. It also im- 
plies the mathematical function approximation for any linear 
and nonlinear systems [3, 5]. 

Gibbs phase rule allows one to determine the number of 
degrees of freedom or variance of a system. This is useful for 
interpreting phase diagrams. 


F=2+C-p, (1) 


where F is the number of degrees of freedom, C is the num- 
ber of chemical components and p is the number of phases in 
the system. The number two is specified because this formu- 
lation assumes that both p and T can be varied [6]. According 
to Eq. (1) in the two-phase region, corresponding to the satu- 
ration curve in the (p— T) diagram, the thermodynamic prop- 
erties of water and steam depend on a single state variable, 
temperature in our case. The range of prediction is valid along 
the entire vapor-liquid saturation line from 273.15 K to the 
critical temperature Tg, i.e., 273.150K < T < 647.096 K [4]. 
The training of neural network is carried out by presenting 
a series of input data and target output values using experi- 
mental database including the entire range of validity where 
temperature is the single input parameter [4, 7-11]. 


Il. ARTIFICIAL NEURAL NETWORKS APPROACH 


Neural networks operate as a black box model requiring 
no detailed information about the system. It is an advanced 
mathematical modeling procedure inspired by biological neu- 
ron systems. The ANN approach seems to be completely suit- 
able to the problems where the relations between variables are 
not linear and complex. In a multi-layer structure (Fig. 1), the 
neurons are grouped into layers, a layer of input neurons, a 
layer of output neurons and one or more hidden layers which 
are made up of many interconnected neurons. 


Hidden 
layer 


Input 
bı layer i 
e e 


Output 
layer 


i=20 
Fig. 1. (Color online) Structure of the ANN. 


The normalization of values is a crucial step in the ANNs. 
The input values to the ANN may differ by several orders 
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of magnitude, which may not reflect the relative importance 
of the inputs determining outlet thermophysical and transport 
properties. To this aim, the input data are normalized within 
the range of [—1, 1] using a mapminmax algorithm, given by 
Eq. (2), to normalize the maximum and minimum values of 
each row. Against the output variables are normalized by log 
function [12]. 


y= (Ymax = Ymin) (x a Lmin) ais Ymin- (2) 


Tmax — Tmin 


After examining a considerable number of differently 
structured neural networks, the adequate ANNs selected in 
this paper have a single hidden layer with 20 neurons and an 
output layer with 10 neurons. The hidden layer has a tansig 
transfer function. The output layer has a purelin transfer func- 
tion. Typical structure of the ANN is shown in Fig. 1. Single 
input variable (T°: temperature) is the input, whereas 10 ther- 
modynamic and transport properties of the density, specific 
heat, dynamic viscosity, thermal conductivity, and speed of 
sound, on saturated liquid and saturated vapor, are considered 
as the output variables. 

Each neuron sums the product of each connection weight 
(wj) from a neuron (j) to the neuron (k) and input (x,;), and 
the additional weight called the bias to get the value of sum 
for the neuron. The i® neuron has a summer that gathers its 
weighted input w;; - x; and the bias b; to form its net input 
P; [3]. 


P= 5 WijXj — bi, (3) 


where w;; denotes the strength of connection from the 7 th in- 
put to the i® neuron, x, is the input vector; b; is the i neu- 
ron bias. An activation function F'(P;), the sigmoid function, 
is used to calculate the neuron output given the set of neu- 
ron inputs. To find suitable ws and biases for each neuron, 
a process training is essential; it is the first step to build an 
ANN. Training means that the weights are corrected to pro- 
duce prespecified (“correct”, known from experiments) target 
values, and the training requires sets of pairs (Xs, Ys) for 
input: the actual input into the network is a vector (Xs), and 
the corresponding target is labeled (Ys) after successful train- 
ing. When correct values of Ys for each vector of Xs from 
the training set are obtained, it is hoped that the network will 
give correct predictions of Y for any new object of X accord- 
ing to the ANN model fundamentals and with the use of more 
data for training the network, better result would be obtained. 
The most utilized training method for multilayered neural net- 
work is called back propagation, where Levenberg-Marquardt 
(LM) is applied which is considered the most efficient algo- 
rithm in terms of speed and memory usage [13]. 

The number of observed data used in the ANN is 377 which 
are divided into three sections: the training set (275 data), test 
set (53 data) and validation set (49 data). Training, test and 
validation subsets of the ANN are obtained as selecting 72% 
of the dataset as training, 14% of the dataset as test and 14% 
of the dataset as validation subsets. 
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The differences between observed and predicted values are 
filtered back through the system and is used to adjust the con- 
nections between the layers, thus performance improves. The 
coefficient of root mean square error (RMSE) is the main cri- 
terion to evaluate the performance of ANN, which is defined 
as follows [3, 5]: 


RMSE = 


Slr 


i 1/2 
X yi- “| (4) 


i=1 


Statistical quality of the ANN for the training, test and 
validation sets is evaluated using the squared correlation 
coefficient R, absolute error AF and average absolute error 
AAE, 


n pat 2 
R=1 diet (yi yi) (5) 
X= (Yi — Yo) 
with 
1% ; 
y= (yi-yi), (6) 
į{=1 
tj A 
M 
1 n 
AAE = — AE, 
> >, (8) 


where y; represents either the i trained, test or validation 
output value and yt is the corresponding target value, with n 
being the number of input vectors. The results are summa- 
rized in Table 1. 

Observed and predicted (dynamic viscosity on saturated 
vapor as example) breakthrough curves shown in Fig. 2 in- 
dicate that the ANN describes the experimental data well. 

Seen to the modeling performances, the general model ob- 
tained from the ANN belonging to all of the properties of 
water and steam was implemented in UDF that compiled and 
hooked in the solver ANSYS-FLUENT. Figure 3 shows the 
incorporation of UDF into the solver of a group of 8 prop- 
erties (density, dynamic viscosity, thermal conductivity and 
speed of sound on saturated liquid and saturated vapor) and 
a group of 2 properties (specific heat on saturated liquid and 
saturated vapor). The UDF written in C++ is a routine which 
can be dynamically linked with the solver FLUENT and pro- 
grammed by the user. 


Ill. CFD CALCULATIONS WITH DEVELOPED UDF 
A. Benchmark case 


The proposed benchmark exercise, experimentally studied 
by Bartolemei and Chanturiya [14], is upward flow of sub- 
cooled water through a heated vertical tube of 615.4m x 2m. 
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Fig. 2. (Color online) Comparison of target and ANN predicted val- 
ues for dynamic viscosity on saturated vapor. 


As shown in Fig. 4, the operating pressure is 4.5 MPa. 
The subcooled water, with the subcooling of 58.2K, en- 
ters from the bottom side and travels upwards through the 
tube. The heat flux applied uniformly on the tube surface 
is 0.57 MW/m?. The inlet mass flux is 900 kg/(m? s). The 
available experimental data include temperatures along the 
tube wall and axis, the bulk liquid temperature, and the cross- 
sectional averaged vapor volume fraction along the tube. The 
wall boiling is in the nucleate boiling regime [15, 16]. 

Parametric study is performed to investigate the effect of 
the developed UDF as summarized in Table 2. 


B. Wall boiling model 


In the nucleate subcooled boiling in a heated tube, the wall 
heat is partially used to form bubbles and the remaining por- 
tion is transferred to the liquid. The heat transfer from the 
wall in the vicinity of a nucleation site occurs during two dis- 
tinct periods: the bubble growth time and the waiting time. 
According to the RPI, the total heat flux from a wall to liquid 
is partitioned into three parts [2, 17] 


wall = de + da =F de, (9) 


where qc is the single-phase convective heat flux, gg is 
quenching heat flux transferred to the liquid phase during the 
waiting time, and dg is the heat flux associated with phase 
change, i.e., evaporation. 


de = he(Ta — Ti)(1 — Ao), (10) 


where h, is the liquid phase heat transfer coefficient; Ty and 
T are the wall and liquid temperature near the wall, respec- 
tively; and A, is the portion covered by nucleation bubble. 


(Te — Ti) Av, (1) 
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TABLE 1. Statistical performance of the trained ANN 


Thermodynamic and transport properties Best linear fit 


Density(liquid) Y = X — 0.03100 
Density(vapor) Y = X — 1.20000 
Specific heat(liquid) Y = X + 0.02500 
Specific heat(vapor) Y = X + 0.01600 
Dynamic viscosity(liquid) Y = X — 0.00800 
Dynamic viscosity(vapor) Y = X — 0.00018 
Thermal conductivity(liquid) Y = X + 0.39000 
Thermal conductivity(vapor) Y = X + 0.00150 
Speed of sound(liquid) Y = X — 0.00970 
Speed of sound(vapor) Y = X — 0.00370 
TABLE 2. Base case values 
Parameters” Value 
Length (m) 2 
Radius (mm) 7.7 
Operating pressure (MPa) 4.5 
Inlet mass flux (kg/(s m?)) 900 
Heat flux (MW/m”) 0.57 
Inlet temperature (K) 473.15 
Inlet subcooling (K) 58.2 


4 Wall, stainless-steel; fluid, H,O 


where kı and A; are heat conductivity and diffusivity in the 
liquid phase and T is the period of bubble detachment. 


{E = VaNapsh fyf, 


where Va is the volume of the bubble based on the bubble 
departure diameter, Ny is the active nucleate site density, py 
is the vapor density, h fy is the latent heat of evaporation, and 
f is the frequency of bubble departure. 

Equations (9)—(12) require closure parameters with empir- 
ical relationship as frequency of bubble departure, bubble de- 
parture diameter, nucleate site density, etc. 


(12) 


C. UDF-RPI model validation 


As the problem formulation is axisymmetric, the domain 
simulated is only a 2D slice with width equal to the tube 
radius. The code manual recommends using a quadrilat- 
eral computational mesh for Eulerian multiphase model, after 
several attempts to find the best computational meshes, we 
adopted grid with 80 uniform radial elements and 1000 uni- 
form axial elements. Finer mesh density in a contact region 
(near the heated wall) to provide a better distribution of local 
flow parameters and to achieve a stable solution. 

To ensure a fully-developed profile of velocity magnitude 
and turbulence quantities at the inlet, the outlet profiles of 
these quantities generated for a simulated flow field without 
boiling (single phase) will be used as inlet information to the 
boiling (multiphase) simulation. 

The Eulerian-RPI method, which allows in FLUENT the 
modeling of multiple separate, yet interacting phases, is em- 
ployed to predict the distribution of the local flow parameters, 


R RMSE AE max (%) AAE (%) 
1 0.1936 0.2395 0.0681 
1 0.3143 0.6419 0.1518 
1 0.3906 3.2169 0.3598 
1 0.2513 2.4516 0.53 

1 0.0881 1.0165 0.2498 
1 0.0019 2.0948 0.1941 
1 0.0592 1.2867 0.2124 
1 0.0505 0.6982 0.2713 
1 0.0192 0.8621 0.1042 
1 0.0071 0.1805 0.0328 


TABLE 3. Input parameters of subcooled nucleate boiling model 


Input Value 
Solver Time Steady 
Type Pressure based 
Velocity formulation Absolute 
2D space Axisymmetric 
Models Energy Active 
Viscous RNG k-e 
Near wall treatment Enhanced 
Multiphase Eulerian 
Drag Ishii 
Lift Moraga 


Turbulent dispersion 
Turbulence interaction 
Heat transfer 
Interfacial area 
Bubble diameter 
Mass transfer 


Lopez-De-Bertodano 
Toshko-Hassan 
Ranz-Marshall 
Ia-Symmetric 
Sauter-Mean 

RPI boiling 


i.e., the vapor volume fraction, the bubble diameter and liq- 
uid temperature. The turbulence phenomena are described 
by a classical Re Normalisation Group (RNG) k-e model in 
combination with enhanced wall treatment for the near-wall 
treatment. Value of yt = 5 is considered reasonable for the 
enhanced wall treatment approach selected. See Table 3 for a 
more detailed list of input parameters and Table 4 for correla- 
tion of the boiling model used in benchmark case [18]. 

The saturation temperature at selected pressure is 
530.55 K, the subcooled flow boiling model water properties 
are assumed to vary with temperature. It is technically part of 
the UDF we developed. That is the aim of this work. 

The governing equations are non linear and several itera- 
tions of loop must be performed before a convergent solution 
is obtained. For the numerical accuracy, the first order up- 
wind scheme is used for spatial discretization of set of gov- 
erning equations. As the numerical scheme employed is cou- 
pled with volume fractions, a low courant number of 10 can 
be used to achieve faster convergence. 


IV. RESULTS AND DISCUSSION 


To validate the proposed UDF in conjunction with RPI 
wall boiling model, comparisons were made among the 
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Fig. 3. (Color online) Developed UDF into CFD solver ANS YS-FLUENT. (a) Group of 8 properties; (b) Group of 2 properties. 


RPI model, the RPI+UDF model, and experimental data 
of Bartolomei and Chanturiya for a heated vertical tube 
at 4.5MPa [14]. The two models were with defining 
solution-dependent material properties as piecewise-linear 
functions of temperature using two data points at 473.15K 
and 543.15 K. Figure 5(a) shows that the proposed UDF-RPI 


axial vapor volume fraction distribution agrees well with ex- 
perimental data [14]. Idem for the liquid temperature profile 
along the tube. Figure 5(b) shows that the proposed UDF- 
RPI is in good agreement with the experimentally measured 
axial liquid temperature distributions. A comparison of the 
predictions with experimental data is demonstrated. 
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TABLE 4. Input parameters of subcooled nucleate boiling model 


Bubble departure diameter 
Tolubinski-Kostanchuk 


Frequency of bubble departure 
Cole 


Subcooled 


water 014.5 mm x2 m 


Tin = Tsat-58.2 K p = 4.54 MPa 


gw=0.57 MW:m-2 


Fig. 4. (Color online) Geometry and flow conditions. 


[= RPI / 
[SRBI FODE PO j 


2 e 
P o 


© 
N 


Vapor volume fraction 
o 
w 


540 
530 PEE aas ia 
x a 
g 520 zo | 
5 ; 
© 510 
oO 
a 
E 500 ee Se 
= =- RPI 
2 
ae H RPIFUDE! 
ea | 
480 |=- Exper, 
470 P ee en ee O es Se ee ee O O BN) 
0.0 0.5 1.0 1.5 2.0 


z/m 


Fig. 5. (Color online) Comparisons of axial vapor volume fraction 
profile (a) and axial liquid temperature (b). 


Parallel to works devoted to axial profiles of vapor vol- 
ume fraction and liquid temperature, the radial distribution 
at various axial locations along the tube was evaluated. Fig- 
ure 6 shows the radial profiles of vapor volume fraction and 
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Fig. 6. (Color online) Profiles of vapor volume fraction (a) and liq- 
uid temperature (b) at 5 locations along the tube. 


liquid temperature at 5 locations along the tube, for the pro- 
posed UDF in conjunction with RPI wall boiling model [19]. 
In Fig. 6(a), boiling is concentrated near the heated wall 
for the flow regime considered here. It is seen that the va- 
por volume fraction is substantially higher adjacent to the 
wall and decreases towards the center. Radial temperature 
profiles at 5 locations are shown in Fig. 6(b). The ther- 
mal non-equilibrium is evident. There are large tempera- 
ture differences between the center and wall of the pipe at 
z = 1.1m, where bubble production starts close to wall while 
the center temperature is below the saturation temperature. 
Non-equilibrium conditions decrease along the pipe and a 
radially uniform temperature distribution is reached around 
z=1.7m. 

As a consequence of the comparison of the present study, 
it is clear that the case of UDF in conjunction with RPI wall 
boiling model, based on artificial neural networks is rea- 
sonable. This demonstrates the strong non-linearity of the 
forms of the thermophysical properties of water and steam of- 
ten requires polynomial interpolations for defining solution- 
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dependent material properties, what makes the difference 
between the two models as compared to the experimental data 
of Bartolomei and Chanturiya, as shown in Fig. 5. 


V. CONCLUSION 


This work is focused on the numerical study of nucle- 
ate subcooled boiling in a heated tube using RPI Wall boil- 
ing model implemented in the solver ANS YS-FLUENT ver- 
sion 14.5.0 in conjunction with the developed UDF based on 
artificial neural networks that calculate the thermodynamic 
and transport properties on saturated liquid and saturated va- 
por of the water assumed to vary with temperature. The val- 
idation of the results has been made with the experimental 
data performed with the baseline case conditions. The re- 
sults demonstrate that the RPI Neural-UDF method is able 
to predict reasonably well the vapor volume fraction distri- 
butions and liquid temperature profile in the heated vertical 
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